\chapter{数值优化}
比如对于某个似然函数$ \mathcal{L}(\bm{\theta}) $，如何找到令$ \mathcal{L} $最大的$ \bm{\theta} $值？前提是我们只有能力根据给定的$ \bm{\theta} $计算$ \mathcal{L}(\bm{\theta}) $。

\section{格点搜索}
假设我们有如下一个似然函数，
\[ \mathcal{L}(\phi) = -\frac{T}{2}\ln 2\pi +\frac{1}{2}\ln (1-\phi^2)-\frac{1}{2}(1-\phi^2)y_1^2-\frac{1}{2}\sum_{i=2}^T(y_t-\phi y_{t-1})^2 \]

其中，$ \phi $是参数，$ y_1,y_2,\cdots,y_T $是已知的观测值。怎样的$ \phi $会使得$ \mathcal{L} $最大？如果我们有5个观测值，$ y_1=0.8,y_2=0.2,y_3=-1.2,y_4=-0.4,y_5=0 $，那么我们可以任意猜测一个$ \phi=0 $，然后计算$ \mathcal{L}(\phi)= -5.73$。依次类推，我们可以在$ (-0.9,0.9) $之间每隔0.1的$ \phi $，分别计算$ \mathcal{L}(\phi) $，可以从中找到一个具有$ \mathcal{L}(\phi) $的$ \phi $。如下图所示，

\begin{figure}[H]\centering
\includegraphics[scale=0.7]{lik_val.png}
\caption{不同$ \phi $下的对数似然值}
\end{figure}

注意：
\begin{itemize}
	\item 格点搜索适用于单峰，如果是多峰时，格子要足够细才能露出峰。
\end{itemize}

\section{黄金分割算法}
\begin{figure}[H]\centering
	\includegraphics[scale=0.7]{golden.pdf}
	\caption{黄金分割算法示意图}
\end{figure}

这个算法是针对单峰函数求最大值。比如在$ [0,1] $之间找函数的最大值，那么可以找两个点$ x_1,x_2 $，如果$ f(x_1)>f(x_2) $则最大值必在$ [0,x_2] $之间，否则最大值在$ [x_1,1] $之间。一旦确定最大值在$ [0,x_2] $之间，则可以类似地再找两个点$ x_1g,x_2g $，继续划分下去。直到想要的精度。

这个算法的问题在于如何确定$ x_1,x_2 $的位置。那么有两个准则：
\begin{itemize}
	\item 最好的办法是要$ [0,x_2] $和$ [x_1,1] $的长度是相同的，这样就不会使得区间长短不一，给我们接下来的搜寻带来不必要的浪费。这个条件要求，
	\[ x_1=1-x_2 \]
	\item 同时，要使得$ x_1$或$x_2$成为下一次区间划分时的一个点，即，
	\[ \frac{1-x_2}{1}=\frac{x_2-x_1}{x_2}\Longrightarrow x_2^2=x_1 \]
\end{itemize}
这两个条件意味着，$ x_2^2+x_2-1=0\Longrightarrow x_2=g=(\sqrt{5}-1)/2=0.618 $。

对于一个区间$ [a,b] $，算法伪码可以书写如下，

\begin{lstlisting}[language = R, mathescape=false]
for i =1,2,3,...
	g = 0.618
	if f(a+(1-g)(b-a)) < f(a+g(b-a)) # f(x1) < f(x2)
		b = a + g(b-a) # x2变成原来的x1
	else
		a = a + (1-g)(b-a) # x1变成原来的x2
	end
end	
\end{lstlisting}

\section{最陡爬坡法}
一个参数可以用格点搜索法，对于多个参数，如果函数连续可微，则可用最陡爬坡法。

假设我们已经有了对参数一个初始猜测，记作$ \bm{\theta}^{(0)} $，然后我们想在一个固定的邻域内找到另外一个$ \bm{\theta}^{(1)} $，它可以使得函数值变得更大。这在数学上可以写成，
\begin{align*}
	\max_{\bm{\theta}^{(1)}}\quad &\mathcal{L}(\bm{\theta}^{(1)})\\	
s.t.\quad &(\bm{\theta}^{(1)}-\bm{\theta}^{(0)})'(\bm{\theta}^{(1)}-\bm{\theta}^{(0)})=k
\end{align*}

对该最优规划可以构造拉格朗日函数如下，
\[ J(\bm{\theta}^{(1)}) = \mathcal{L}(\bm{\theta}^{(1)})+\lambda [k-(\bm{\theta}^{(1)}-\bm{\theta}^{(0)})'(\bm{\theta}^{(1)}-\bm{\theta}^{(0)})]\]

对上式进求导并令其为0，有，
\[\left. \frac{\partial \mathcal{L}(\bm{\theta})}{\partial\bm{\theta}}\right|_{\bm{\theta}=\bm{\theta}^{(1)}}-2\lambda (\bm{\theta}^{(1)}-\bm{\theta}^{(0)})=0 \]

上式可以进一步写成，
\[ \bm{\theta}^{(1)}-\bm{\theta}^{(0)}=\frac{1}{2\lambda}\cdot \left. \frac{\partial \mathcal{L}(\bm{\theta})}{\partial\bm{\theta}}\right|_{\bm{\theta}^{(1)}}\]

而$ \left. \frac{\partial \mathcal{L}(\bm{\theta})}{\partial\bm{\theta}}\right|_{\bm{\theta}^{(1)}} $实际上是我们所说的梯度，可以把它记成$ g(\bm{\theta}^{(1)}) $。这样，最优的$ \bm{\theta}^{(1)} $就可以通过下式找到，
\begin{equation}\label{opt_eq_climb}
\bm{\theta}^{(1)}-\bm{\theta}^{(0)}=\frac{1}{2\lambda}\cdot g(\bm{\theta}^{(1)})	
\end{equation}


注意到：
\begin{itemize}
	\item 仅允许$ \bm{\theta} $变动一个固定值，那么最优的变动为$ 1/2\lambda $乘以梯度向量$ g(\bm{\theta}^{(1)}) $。$ \lambda $意味着某个正数\footnote{利用一阶条件和约束条件是可以计算出$ \lambda $的，但实践中我们从不会去计算，只是设定一个小小的正数即可。}。
	\item 如果只变动很小的一步，那么$ g(\bm{\theta}^{(1)}) $ 与$ g(\bm{\theta}^{(0)}) $将非常接近，实践中通常用$ g(\bm{\theta}^{(0)}) $替代$ g(\bm{\theta}^{(1)}) $进行计算。
\end{itemize}

\subsection{一个数值例子}
如果只有2个参数，令似然函数为，
\[ \mathcal{L} = -1.5\theta_1^2-2\theta_2^2 \]
上式的梯度向量为，
\[\begin{bmatrix}
	\frac{\partial \mathcal{L}(\bm{\theta})}{\partial \theta_1}=-3\theta_1\\
		\frac{\partial \mathcal{L}(\bm{\theta})}{\partial \theta_2}=-4\theta_2
\end{bmatrix}
\]

初始点猜测为$ \bm{\theta}^{(0)}=(-1,1)' $，那么此时梯度向量的值为，
\[\begin{bmatrix}
	\left.\frac{\partial \mathcal{L}(\bm{\theta})}{\partial \theta_1}\right|_{\bm{\theta}=\bm{\theta}^{(0)}}=-3\theta_1=3\\
	\left.\frac{\partial \mathcal{L}(\bm{\theta})}{\partial \theta_2}\right|_{\bm{\theta}=\bm{\theta}^{(0)}}=-4\theta_2=-4
\end{bmatrix}
\]
那么，当$ k=1 $时，就可以反算出$ \lambda $，从而依据\eqref{opt_eq_climb}式，有，
\begin{align*}
	\theta_1^{(1)}-	\theta_1^{(0)} =\frac{3}{5}\\
		\theta_2^{(1)}-	\theta_2^{(0)} =-\frac{4}{5}\\
\end{align*}
然后重复这一过程，直到想要的精度。一个双变量的最陡爬坡示意图如下所示。当然，该法找到的也是局部最大值，这取决于初值的设定。
\begin{figure}[H]
\includegraphics[scale=0.7]{climb.png}
\end{figure}
也可以参考\lstinline|contour.mw|文件中的直观图。

\section{牛顿-拉夫森方法(Newton-Raphson)}
最陡爬坡需要大量迭代，NR方法收敛会更快。但它有两个条件：
\begin{itemize}
	\item 对数似然函数二阶导数存在。
	\item 对数似然函数是凹的，即海塞矩阵负定。
\end{itemize}

先看NR法的原理，对似然函数$ \mathcal{L}(\bm{\theta}) $在初值$ \bm{\theta}^{(0)} $处二阶泰勒展开，有，
\begin{equation}\label{opt_eq_nr}
\mathcal{L}(\bm{\theta}) \approx \mathcal{L}(\bm{\theta}^{(0)})+g'(\bm{\theta}^{(0)})(\bm{\theta}-\bm{\theta}^{(0)})-\frac{1}{2}(\bm{\theta}-\bm{\theta}^{(0)})'H(\bm{\theta}^{(0)})(\bm{\theta}-\bm{\theta}^{(0)})	
\end{equation}


其中，
\[\left. \underbrace{H(\bm{\theta}^{(0)})}_{a\times a} = -\frac{\partial^2 \mathcal{L}(\bm{\theta})}{\partial\bm{\theta}\partial\bm{\theta}'}\right|_{\bm{\theta}=\bm{\theta}^{(0)}}\]

对\eqref{opt_eq_nr}式针对$ \bm{\theta} $最优化，即一阶求导，并令之为0，有，
\[ g(\bm{\theta}^{(0)})-H(\bm{\theta}^{(0)})(\bm{\theta}-\bm{\theta}^{(0)})=0 \]
这意味着新的估计$ \bm{\theta}^{(1)} $为，
\[\bm{\theta}^{(1)}-\bm{\theta}^{(0)} =  [H(\bm{\theta}^{(0)})]^{-1}g(\bm{\theta}^{(0)}) \]
同时，注意到我们最优化时没有设置步长，如果与最陡爬坡一样设置了步长的话，新的估计应该是，
\[\bm{\theta}^{(1)} = \bm{\theta}^{(0)}+ s\cdot [H(\bm{\theta}^{(0)})]^{-1}\cdot g(\bm{\theta}^{(0)}) \]

但要注意，如果似然函数不是凹的，那么NR法不如最陡爬坡。

此外，由于二阶导数$ [H(\bm{\theta}^{(0)})]^{-1} $的计算非常费时，后续提出了许多修正的NR方法。